/*
 * Copyright 2013 Google Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package com.aplink.generic.google.maps;

import java.util.ArrayList;
import java.util.List;

import com.google.android.gms.maps.model.LatLng;

public class PolyUtil {

	private static final double DEFAULT_TOLERANCE = 0.1; // meters.

	/**
	 * Computes whether the given point lies inside the specified polygon. The
	 * polygon is always cosidered closed, regardless of whether the last point
	 * equals the first or not. Inside is defined as not containing the South
	 * Pole -- the South Pole is always outside. The polygon is formed of great
	 * circle segments if geodesic is true, and of rhumb (loxodromic) segments
	 * otherwise.
	 */
	public static boolean containsLocation(final LatLng point,
			final List<LatLng> polygon, final boolean geodesic) {
		final int size = polygon.size();
		if (size == 0) {
			return false;
		}
		final double lat3 = Math.toRadians(point.latitude);
		final double lng3 = Math.toRadians(point.longitude);
		final LatLng prev = polygon.get(size - 1);
		double lat1 = Math.toRadians(prev.latitude);
		double lng1 = Math.toRadians(prev.longitude);
		int nIntersect = 0;
		for (final LatLng point2 : polygon) {
			final double dLng3 = MathUtil.wrap(lng3 - lng1, -Math.PI, Math.PI);
			// Special case: point equal to vertex is inside.
			if ((lat3 == lat1) && (dLng3 == 0)) {
				return true;
			}
			final double lat2 = Math.toRadians(point2.latitude);
			final double lng2 = Math.toRadians(point2.longitude);
			// Offset longitudes by -lng1.
			if (PolyUtil.intersects(lat1, lat2,
					MathUtil.wrap(lng2 - lng1, -Math.PI, Math.PI), lat3, dLng3,
					geodesic)) {
				++nIntersect;
			}
			lat1 = lat2;
			lng1 = lng2;
		}
		return (nIntersect & 1) != 0;
	}

	/**
	 * Decodes an encoded path string into a sequence of LatLngs.
	 */
	public static List<LatLng> decode(final String encodedPath) {
		final int len = encodedPath.length();

		// For speed we preallocate to an upper bound on the final length, then
		// truncate the array before returning.
		final List<LatLng> path = new ArrayList<LatLng>();
		int index = 0;
		int lat = 0;
		int lng = 0;

		while (index < len) {
			int result = 1;
			int shift = 0;
			int b;
			do {
				b = encodedPath.charAt(index++) - 63 - 1;
				result += b << shift;
				shift += 5;
			} while (b >= 0x1f);
			lat += (result & 1) != 0 ? ~(result >> 1) : (result >> 1);

				result = 1;
				shift = 0;
				do {
					b = encodedPath.charAt(index++) - 63 - 1;
					result += b << shift;
					shift += 5;
				} while (b >= 0x1f);
				lng += (result & 1) != 0 ? ~(result >> 1) : (result >> 1);

					path.add(new LatLng(lat * 1e-5, lng * 1e-5));
		}

		return path;
	}

	/**
	 * Encodes a sequence of LatLngs into an encoded path string.
	 */
	public static String encode(final List<LatLng> path) {
		long lastLat = 0;
		long lastLng = 0;

		final StringBuffer result = new StringBuffer();

		for (final LatLng point : path) {
			final long lat = Math.round(point.latitude * 1e5);
			final long lng = Math.round(point.longitude * 1e5);

			final long dLat = lat - lastLat;
			final long dLng = lng - lastLng;

			PolyUtil.encode(dLat, result);
			PolyUtil.encode(dLng, result);

			lastLat = lat;
			lastLng = lng;
		}
		return result.toString();
	}

	private static void encode(long v, final StringBuffer result) {
		v = v < 0 ? ~(v << 1) : v << 1;
		while (v >= 0x20) {
			result.append(Character.toChars((int) ((0x20 | (v & 0x1f)) + 63)));
			v >>= 5;
		}
		result.append(Character.toChars((int) (v + 63)));
	}

	/**
	 * Computes whether the vertical segment (lat3, lng3) to South Pole
	 * intersects the segment (lat1, lng1) to (lat2, lng2). Longitudes are
	 * offset by -lng1; the implicit lng1 becomes 0.
	 */
	private static boolean intersects(final double lat1, final double lat2,
			final double lng2, final double lat3, final double lng3,
			final boolean geodesic) {
		// Both ends on the same side of lng3.
		if (((lng3 >= 0) && (lng3 >= lng2)) || ((lng3 < 0) && (lng3 < lng2))) {
			return false;
		}
		// Point is South Pole.
		if (lat3 <= (-Math.PI / 2)) {
			return false;
		}
		// Any segment end is a pole.
		if ((lat1 <= (-Math.PI / 2)) || (lat2 <= (-Math.PI / 2))
				|| (lat1 >= (Math.PI / 2)) || (lat2 >= (Math.PI / 2))) {
			return false;
		}
		if (lng2 <= -Math.PI) {
			return false;
		}
		final double linearLat = ((lat1 * (lng2 - lng3)) + (lat2 * lng3))
				/ lng2;
		// Northern hemisphere and point under lat-lng line.
		if ((lat1 >= 0) && (lat2 >= 0) && (lat3 < linearLat)) {
			return false;
		}
		// Southern hemisphere and point above lat-lng line.
		if ((lat1 <= 0) && (lat2 <= 0) && (lat3 >= linearLat)) {
			return true;
		}
		// North Pole.
		if (lat3 >= (Math.PI / 2)) {
			return true;
		}
		// Compare lat3 with latitude on the GC/Rhumb segment corresponding to
		// lng3.
		// Compare through a strictly-increasing function (tan() or mercator())
		// as convenient.
		return geodesic ? Math.tan(lat3) >= PolyUtil.tanLatGC(lat1, lat2, lng2,
				lng3) : MathUtil.mercator(lat3) >= PolyUtil.mercatorLatRhumb(
						lat1, lat2, lng2, lng3);
	}

	/**
	 * Same as {@link #isLocationOnEdge(LatLng, List, boolean, double)} with a
	 * default tolerance of 0.1 meters.
	 */
	public static boolean isLocationOnEdge(final LatLng point,
			final List<LatLng> polygon, final boolean geodesic) {
		return PolyUtil.isLocationOnEdge(point, polygon, geodesic,
				PolyUtil.DEFAULT_TOLERANCE);
	}

	/**
	 * Computes whether the given point lies on or near the edge of a polygon,
	 * within a specified tolerance in meters. The polygon edge is composed of
	 * great circle segments if geodesic is true, and of Rhumb segments
	 * otherwise. The polygon edge is implicitly closed -- the closing segment
	 * between the first point and the last point is included.
	 */
	public static boolean isLocationOnEdge(final LatLng point,
			final List<LatLng> polygon, final boolean geodesic,
			final double tolerance) {
		return PolyUtil.isLocationOnEdgeOrPath(point, polygon, true, geodesic,
				tolerance);
	}

	private static boolean isLocationOnEdgeOrPath(final LatLng point,
			final List<LatLng> poly, final boolean closed,
			final boolean geodesic, final double toleranceEarth) {
		final int size = poly.size();
		if (size == 0) {
			return false;
		}
		final double tolerance = toleranceEarth / MathUtil.EARTH_RADIUS;
		final double havTolerance = MathUtil.hav(tolerance);
		final double lat3 = Math.toRadians(point.latitude);
		final double lng3 = Math.toRadians(point.longitude);
		final LatLng prev = poly.get(closed ? size - 1 : 0);
		double lat1 = Math.toRadians(prev.latitude);
		double lng1 = Math.toRadians(prev.longitude);
		if (geodesic) {
			for (final LatLng point2 : poly) {
				final double lat2 = Math.toRadians(point2.latitude);
				final double lng2 = Math.toRadians(point2.longitude);
				if (PolyUtil.isOnSegmentGC(lat1, lng1, lat2, lng2, lat3, lng3,
						havTolerance)) {
					return true;
				}
				lat1 = lat2;
				lng1 = lng2;
			}
		} else {
			// We project the points to mercator space, where the Rhumb segment
			// is a straight line,
			// and compute the geodesic distance between point3 and the closest
			// point on the
			// segment. This method is an approximation, because it uses
			// "closest" in mercator
			// space which is not "closest" on the sphere -- but the error is
			// small because
			// "tolerance" is small.
			final double minAcceptable = lat3 - tolerance;
			final double maxAcceptable = lat3 + tolerance;
			double y1 = MathUtil.mercator(lat1);
			final double y3 = MathUtil.mercator(lat3);
			final double[] xTry = new double[3];
			for (final LatLng point2 : poly) {
				final double lat2 = Math.toRadians(point2.latitude);
				final double y2 = MathUtil.mercator(lat2);
				final double lng2 = Math.toRadians(point2.longitude);
				if ((Math.max(lat1, lat2) >= minAcceptable)
						&& (Math.min(lat1, lat2) <= maxAcceptable)) {
					// We offset longitudes by -lng1; the implicit x1 is 0.
					final double x2 = MathUtil.wrap(lng2 - lng1, -Math.PI,
							Math.PI);
					final double x3Base = MathUtil.wrap(lng3 - lng1, -Math.PI,
							Math.PI);
					xTry[0] = x3Base;
					// Also explore wrapping of x3Base around the world in both
					// directions.
					xTry[1] = x3Base + (2 * Math.PI);
					xTry[2] = x3Base - (2 * Math.PI);
					for (final double x3 : xTry) {
						final double dy = y2 - y1;
						final double len2 = (x2 * x2) + (dy * dy);
						final double t = len2 <= 0 ? 0 : MathUtil.clamp(
								((x3 * x2) + ((y3 - y1) * dy)) / len2, 0, 1);
						final double xClosest = t * x2;
						final double yClosest = y1 + (t * dy);
						final double latClosest = MathUtil
								.inverseMercator(yClosest);
						final double havDist = MathUtil.havDistance(lat3,
								latClosest, x3 - xClosest);
						if (havDist < havTolerance) {
							return true;
						}
					}
				}
				lat1 = lat2;
				lng1 = lng2;
				y1 = y2;
			}
		}
		return false;
	}

	/**
	 * Same as {@link #isLocationOnPath(LatLng, List, boolean, double)}
	 *
	 * with a default tolerance of 0.1 meters.
	 */
	public static boolean isLocationOnPath(final LatLng point,
			final List<LatLng> polyline, final boolean geodesic) {
		return PolyUtil.isLocationOnPath(point, polyline, geodesic,
				PolyUtil.DEFAULT_TOLERANCE);
	}

	/**
	 * Computes whether the given point lies on or near a polyline, within a
	 * specified tolerance in meters. The polyline is composed of great circle
	 * segments if geodesic is true, and of Rhumb segments otherwise. The
	 * polyline is not closed -- the closing segment between the first point and
	 * the last point is not included.
	 */
	public static boolean isLocationOnPath(final LatLng point,
			final List<LatLng> polyline, final boolean geodesic,
			final double tolerance) {
		return PolyUtil.isLocationOnEdgeOrPath(point, polyline, false,
				geodesic, tolerance);
	}

	private static boolean isOnSegmentGC(final double lat1, final double lng1,
			final double lat2, final double lng2, final double lat3,
			final double lng3, final double havTolerance) {
		final double havDist13 = MathUtil.havDistance(lat1, lat3, lng1 - lng3);
		if (havDist13 <= havTolerance) {
			return true;
		}
		final double havDist23 = MathUtil.havDistance(lat2, lat3, lng2 - lng3);
		if (havDist23 <= havTolerance) {
			return true;
		}
		final double sinBearing = PolyUtil.sinDeltaBearing(lat1, lng1, lat2,
				lng2, lat3, lng3);
		final double sinDist13 = MathUtil.sinFromHav(havDist13);
		final double havCrossTrack = MathUtil
				.havFromSin(sinDist13 * sinBearing);
		if (havCrossTrack > havTolerance) {
			return false;
		}
		final double havDist12 = MathUtil.havDistance(lat1, lat2, lng1 - lng2);
		final double term = havDist12 + (havCrossTrack * (1 - (2 * havDist12)));
		if ((havDist13 > term) || (havDist23 > term)) {
			return false;
		}
		if (havDist12 < 0.74) {
			return true;
		}
		final double cosCrossTrack = 1 - (2 * havCrossTrack);
		final double havAlongTrack13 = (havDist13 - havCrossTrack)
				/ cosCrossTrack;
		final double havAlongTrack23 = (havDist23 - havCrossTrack)
				/ cosCrossTrack;
		final double sinSumAlongTrack = MathUtil.sinSumFromHav(havAlongTrack13,
				havAlongTrack23);
		return sinSumAlongTrack > 0; // Compare with half-circle == PI using
		// sign of sin().
	}

	/**
	 * Returns mercator(latitude-at-lng3) on the Rhumb line (lat1, lng1) to
	 * (lat2, lng2). lng1==0.
	 */
	private static double mercatorLatRhumb(final double lat1,
			final double lat2, final double lng2, final double lng3) {
		return ((MathUtil.mercator(lat1) * (lng2 - lng3)) + (MathUtil
				.mercator(lat2) * lng3)) / lng2;
	}

	/**
	 * Returns sin(initial bearing from (lat1,lng1) to (lat3,lng3) minus initial
	 * bearing from (lat1, lng1) to (lat2,lng2)).
	 */
	private static double sinDeltaBearing(final double lat1, final double lng1,
			final double lat2, final double lng2, final double lat3,
			final double lng3) {
		final double sinLat1 = Math.sin(lat1);
		final double cosLat2 = Math.cos(lat2);
		final double cosLat3 = Math.cos(lat3);
		final double lat31 = lat3 - lat1;
		final double lng31 = lng3 - lng1;
		final double lat21 = lat2 - lat1;
		final double lng21 = lng2 - lng1;
		final double a = Math.sin(lng31) * cosLat3;
		final double c = Math.sin(lng21) * cosLat2;
		final double b = Math.sin(lat31)
				+ (2 * sinLat1 * cosLat3 * MathUtil.hav(lng31));
		final double d = Math.sin(lat21)
				+ (2 * sinLat1 * cosLat2 * MathUtil.hav(lng21));
		final double denom = ((a * a) + (b * b)) * ((c * c) + (d * d));
		return denom <= 0 ? 1 : ((a * d) - (b * c)) / Math.sqrt(denom);
	}

	/**
	 * Returns tan(latitude-at-lng3) on the great circle (lat1, lng1) to (lat2,
	 * lng2). lng1==0. See http://williams.best.vwh.net/avform.htm .
	 */
	private static double tanLatGC(final double lat1, final double lat2,
			final double lng2, final double lng3) {
		return ((Math.tan(lat1) * Math.sin(lng2 - lng3)) + (Math.tan(lat2) * Math
				.sin(lng3))) / Math.sin(lng2);
	}

	private PolyUtil() {
	}
}
